#pragma once

/**
 * Wrapper around cublasSger()
 * Performs B = alpha * x * y' + A
 *
 * @param alpha float to scale x * y'
 * @param X a length m vector to be scaled and multiplied
 * @param Y a length n vector to be scaled and multiplied
 * @param A a m x n matrix to be added to
 * @param B a m x n matrix to store the result of the computation
 */
void sger
(
	float alpha,
	Vector X,
	Vector Y,
	Matrix A,
	Matrix B
)
{
	bool is_col_A;
	int m, n, lda, size_a;

	is_col_A = (A.stride == 1);

	m = A.length;
	n = A.sub.length;

	if(!is_col_A){
		//warn user
		std::cerr << "sger: WARNING: Passing in a row major matrix has incurred a copy operation" << std::endl;
		int size_A = n * m * sizeof(float);
		float *A_data;
		cudaMalloc(&A_data, size_A);
		SubMatrix A_sub(n, m, A_data, 0);
		Matrix A_copy(m, 1, A_sub);
		copy(A, A_copy);
		A = A_copy; //huzzah for args passed by copy
		is_col_A = 1;
	}

	if(m != X.length || n != Y.length){
		std::cerr << "sger: ERROR length mismatch in X or Y." <<std::endl;
	}

	lda = A.sub.stride;
	size_a = m * lda;
	//copy A into the result buffer
	cudaMemcpy(B.sub.data, A.sub.data, size_a, cudaMemcpyDeviceToDevice);

	cublasSger(m, n, alpha, X.data, X.stride, Y.data, Y.stride, B.sub.data, lda);
}

/**another wrapper that uses stored sequences instead of 
 * uniform nested sequences for vectors
 */

void sger
(
 	float alpha,
	stored_sequence<float> X,
	stored_sequence<float> Y,
	Matrix A,
	Matrix B
)
{
	Vector x(X.length, 1, X.data);
	Vector y(Y.length, 1, Y.data);
	sger(alpha, x, y, A, B);
}

